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Abstract 

Numerical schemes derived from gas-kinetic theory can be applied to 
simulations in the hydrodynamics limit, in laminar and also turbulent 
regimes. In the latter case, the underlying Boltzmann equation describes 
a distribution of eddies, in line with the concept of eddy viscosity devel- 
oped by Lord Kelvin and Osborne Reynolds at the end of the nineteenth 
century. These schemes are physically more consistent than schemes de- 
rived from the Navier-Stokes equations, which invariably assume infinite 
collisions between gas particles (or interactions between eddies) in the 
calculation of advective fluxes. In fact, in continuum regime too, the lo- 
cal Knudsen number can exceed the value 0.001 in shock layers, where 
gas-kinetic schemes outperform Navier-Stokes schemes, as is well known. 

Simulation of turbulent flows benefit from the application of gas-kinetic 
schemes, as the turbulent Knudsen number (the ratio between the eddies' 
mean free path and the mean flow scale) can locally reach values well 
in excess of 0.001, not only in shock layers. A further advantage of gas- 
kinetic schemes is that the fluxes are accurate to r 2 , for instance in the 
scheme developed in [19] for the finite-volume discretization. In laminar 
flow, this provides a better resolution of shocks and vortexes, whereas in 
turbulent flows, high-order fluxes allow for a better resolution of secondary 
flows in a manner comparable to higher-order turbulence models for the 
Navier-Stokes schemes. 

This study has investigated a few cases of shock - boundary layer inter- 
action comparing a gas-kinetic scheme and a Navier-Stokes one, both with 
a standard k — uj turbulence model. Whereas the results obtained from 
the Navier-Stokes scheme are affected by the limitations of eddy viscosity 
two-equation models, the gas-kinetic scheme has performed much better 
without making any further assumption on the turbulent structures. 

1 Introduction 

Numerical schemes based on the Navier-Stokes (NS) equations have benefited 
for many years from the valuable work of mathematicians and engineers. Aero- 
dynamicists dispose of accurate and fast simulation tools, which can be applied 



to complex geometries and challenging flow conditions, providing physically 
consistent results in many cases. 

However, limitations still apply and appear difficult to overcome. They 
concern the numerical model - dependency of the results on numerical scheme 
and mesh, steeply increasing computational cost with the order of the scheme, 
reduction of accuracy at discontinuities - or the physical model - modelling 
the effect of unresolved turbulence on resolved flow in conditions far from local 
turbulent equilibrium, simulation of rarefied flow. The prediction of turbulent, 
hypersonic flow is a particularly challenging example. 

Numerical schemes based on gas-kinetic theory, or Gas-Kinetic Schemes 
(GKS), might help in overcoming these limitations. Firstly, because the Boltz- 
mann equation is a physically more accurate representation of fluid mechanics 
than the Navier-Stokes equations and secondly, because gas-kinetic can deal 
much better with the discontinuities that invariably at cells interface in most 
numerical approaches, since Godunov scheme ( (jj). 

Navier-stokes schemes split transport and collisions, i.e. advective and vis- 
cous fluxes. As such, advective fluxes are generated by the solution of the Rie- 
mann problem which does not consider the effect of particle collision. The effect 
of collisions is added a posteriori by the viscous fluxes, calculated independently. 
As long as the molecular relaxation time is much smaller than the mean time 
scale, superposition of fluxes does not affect the physical consistence. However, 
wherever this assumption is not true, as in the case of shock layers or of rarefied 
flows, the physical consistence of Navier-Stokes schemes becomes questionable. 
Despite the fact that Navier-Stokes schemes can predict shocks satisfactorily for 
many industrial applications, the prediction treats the shock merely as a discon- 
tinuity. Moreover, special treatment is often necessary to maintain the stability 
of the solution in presence of shocks - and this contributes to the dependence of 
the solution from the chosen method. 

The modelling of turbulence can also benefit from the use of GKS. It is rec- 
ognized since the publications of Lord Kelvin and Osborne Reynolds' works on 
turbulence, (refer to [3|[4j and references therein) that the Boltzmann equation, 
also in the simpler form of the BGK model, can be used to describe not only 
the flow as distribution of particles but also the turbulent flow as distribution 
of eddies. Moreover, it is also recognized ( [4]) that the projection onto the 
physical space (x,t) of the BGK model generates a higher-order (in r t ) turbu- 
lent stress tensor. In particular a third order Chapman-Enskog expansion for 
/ generates a turbulent stress tensor of the second order, i.e. comparable to 
non-linear turbulent stress models. These models are known to provide more 
accurate values than linear models ( [l3|[l7] ). 

This study focuses on the GKS developed by Xu ( [19]) which has also 



been investigated by other researchers ( [121 16 ) and has provided very good 



results in a number of cases, ranging from viscous-dominated, subsonic flows to 
hypersonics. It provides fluxes accurate to r 2 and, consequently, a second-order 
turbulent stress tensor. 

This GKS has been implemented into an existing solver which uses a stan- 
dard k — uj ( [18]) turbulence model. The flow cases investigated are popu- 
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lar aerodynamic benchmarks - all in the continuum regime - characterised by 
strong shock - boundary layer interaction - that is a flow condition where most 
two-equation turbulence models fail to predict shock position and extension of 
separated flow accurately. The equation for the turbulent quantities k and u> 
are advanced in a segregate way. 

This paper presents a brief description of the GKS for laminar and turbulent 
flow, followed by the numerical experiments and conclusions. 



2 Gas-kinetic scheme 

A few gas-kinetic schemes for the solution of the Euler and the Navier-Stokes 
equations have been proposed in the 1990s ( [5|[Tl][T2j[l9j[2l] and references 
therein) as an alternative to the most popular schemes, which normally assume 
continuity of the flow or solve a Riemann problem at cells interfaces. 

The main idea behind these schemes is to consider a discontinuous state 
across interfaces, re-construct the equilibrium and non-equilibrium distribution 
functions based on the macroscopic flow variables and calculate the evolution of 
the distribution functions during a time step At integrating the BGK-Boltzmann 
equation. The macroscopic flow quantities arc then recovered taking moments 
of the solution distribution function. 

We use the macroscopic variables p, U = [ui u-i U3] 71 , and E to describe 
density, velocity and total energy of a gas. Instead of using the well-known 
Navier-Stokes equations, we write the BGK model following [l]: 

df df (g - f) 

-d-t +u ^r~r- (1) 

where the summation convention holds, / is the gas distribution function, g is 
the equilibrium state, a Maxwellian distribution, approached by / and r is the 
particles collision time, which is related to the molecular viscosity and heath 



coefficients of the gas. Although not explicitly indicated, it is assumed in 19 
that the collision time can also include the effects of turbulence, beside those of 
molecular viscosity and numerical dissipation. The variable £ is related to the 
additional degrees of freedom of the gas molecules. £ has K degrees of freedom, 
where: 

where 7 is the specific heat ratio. The equilibrium distribution is: 

g = P ^e-^'^ + e) (3) 

where A = m is the molecular mass, k is the Boltzmann constant, and T is 
temperature. The relation between macroscopic variables and gas distribution 
function is: 
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where ip is: 




note that dS = du\ du2 du 3 £ 1 d£. Conservation of mass, momentum and 
energy during particle collision is expressed by: 

J(g-f)^dZ = (6) 
The BGK equation [T] has an analytical solution: 



1 f l 

f(x, y, z, t, u, v, w, f ) = - / g(x', y', z', t, u, v, w, C)e" (t ~' )/r dt'+e~ t/r f (x-ut, y-vt, z-wt) 



(7) 

where fo is the initial gas distribution function, x' — x — u(t — if), y' = y — 
v(t — tf), z' = z — w(t — if). The kernel of the GKS consists in expressing the 
distribution function / at cells or volumes interfaces in order to assess the fluxes 
as functions of /. For instance the flux in direction i at the interface between 
cells n and n + 1 can be expressed as a first moment of /: 

Ff +1/2 = [ < +1/ V l+1/2 f(x n+1 / 2 )d~dt (8) 



(i 



The distribution functions fo and g in the [7] must be consistent with the macro- 
scopic variables and their gradients. An important assumption in the deriva- 
tion of the GKS is that whereas equilibrium distributions are Maxwellians, the 
non-equilibrium distribution are expressed as Taylor expansion of Maxwellian 
distribution. Assuming an interface normal to direction 1 located at x\ = 0, 
the initial equilibrium distribution is expressed as: 

J .9o (l + a\x t -At), x l <Q . . 

g \ g Q (1 + aV Xi -At), Xl >0 w 

where go is a Maxwellian derived from a state [po pU^i PqEq], which is an 
average state between left and right, obtained in a non-trivial averaging process, 
which fulfils the BGK model and the conservation laws ( [l9]). The initial 
distribution f can be expressed as: 

f g l (l + a\ Xi ) -T(a\ui+A l ) , xi < . . 

J0 \ g r (l + a r iXi )-T(a r iUi +A r ), x x >0 K n 
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where g l and g r are Maxwellian distribution on both sides of the interface, which 
are indicated as left and right. The choice of the terms used in the expansion 
[l0| is critical for the type of GKS. The terms proportional to r represent the 
non-equilibrium parts in the Chapman-Enskog expansion ( [2 ), whereas the 
expansion in the spatial directions Xi is directly related to the formal accuracy 
of the resulting scheme. Moreover, one can have a directional splitting scheme 
by simply expanding in the direction normal to the interface, or a truly multi- 
dimensional scheme by considering the derivatives in all directions. 
Each of the coefficients in the [9| and |l0| is expanded as: 

o>% = a-ii + a-i2 u + a i3 v + ai/jW + cii 5 (u 2 + v 2 + w 2 + £ 2 ) (11) 

All the components of the coefficients arc determined from compatibility 
relations with the macroscopic variables and the [6| The details can be found 
in |19| . The determination of all coefficients involves the solution of numerous 
(depending on the dimensions) linear systems and the evaluation of the erfc 
function, which contribute to the computational cost. Inserting the [9] and [10] 
into the[7J we obtain /: 

/ = (l - e~' /T ) ,9o + (-T + re"*/ T + * e^'A (ti oj 1 + h r of) Ui g + (t - r + re^'/j Ag + 
+ e~ t/T (h l g l + h r g r - (t + r) ( Ul a l t h'g l + Ul a\h T g r )) - Te~ t/T (A l h l g l + A r h r g r ) (12) 

where h — H(U), h r = 1 — H(u) and all coefficients are intended as series 



expansions in the form of 11 Advective and viscous fluxes cannot be clearly 
separated in the[l2] this is a consequence of the fact that transport and collision 
of particles / eddies are considered simultaneously. Like in other gas-kinetic 
scheme the resulting fluxes appear as a series expansion in r, in this case accurate 
to the second order. Zero order terms provide the advective fluxes corresponding 
to the average state <?0j first order terms include the contribution to the viscous 
fluxes of the average state go plus a correction to the advective fluxes, whereas 
second-order terms contain correction to the viscous fluxes and represent the 
real higher order contribution. 

The relaxation time r is set as a function of the molecular viscosity of the gas 
plus an additional term proportional to the pressure jump across the interface. 

r=^ + \^4At (13) 

V \p r +p l \ 

A known drawback of the BGK model is that it implies a unity Prandtl number; 
the heath flux must therefore be corrected for realistic gas / fluids ( [19]). 



3 Turbulence modelling 

The present study is based on a simple modelling technique: the Reynolds 
approach, which resolves explicitly the mean flow and models the effects of 
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all turbulent length scales - often referred to as RANS in its implementation 
with the Navier-Stokes equations. Turbulent quantities are modelled according 
to k — uj model ( [18]) two-equation models, which is a popular and accepted 
representative of this class of models. 

OpK dpujK d (. s9K\ „ ., 

— P — /3 puK + — — I (p, + <r Ih)-k— ) (14) 



Ot dxj dxj \ dxj 

Opui dpujUj ui 2 ( , , Qui \ 

~W + J k~- "">K P + ^ + Ox- + ^ Ox- ) (15) 

» t =r p -- (i6) 

UJ 

where P is turbulence production term: 



0~Xj 



Tij is the turbulent stresses tensor: 



s n \ (";'; • n : : ) (i9) 



and Sij is the strain rate: 



2 \Oxj Oxi 
A typical choice for the parameters is: 

P = —, /8* = — , 7=-, 7* = 1, tr=-, <r* = - (20) 
H 40' H 100' ' 9' ' 2' 2 v ' 

The implementation of a turbulence model into a gas-kinetic scheme might seem 

straight-forward and practical steps taken to include the effects of turbulence 

into a GKS-based computation are really simple. Following [19] we can simply 

re-write the BGK model [l] replacing the molecular relaxation time r with a 

relaxation time r* which considers both molecular and turbulent phenomena. 

°l +Ul K = (*zJ) (21) 

Ot 1 OXi T* 

where 

T* =T + T t (22) 

where r t is the turbulent relaxation time, which can be expressed as a linear 
function of turbulent viscosity: 

r t = pk =^ (23) 
ojp p 
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The [T] can be rewritten as: 



In the practical calculations the |13[ used in laminar flow, becomes: 

r* = ^ + ^4At (25) 
p \p r +P l \ 

The turbulent relaxation time can also be expressed as a non-linear function of 
the macroscopic variables as suggested by Chen (adapted from [3]): 

T* = ^ 1/2 (26) 

where rj = S/u is the ratio between unresolved and resolved turbulence time 
scales, where S is a scalar representing local velocity gradient and r is the 
molecular relaxation time To = [ijp- The final expression for the relaxation 
time including numerical dissipation is: 



r* = » + pk/ \ /2 +cf^M (27) 
P p(l + rj 2 ) 1/2 lP r +P l \ 

where the coefficient C is determined heuristically (on average C has been fixed 
at around 0.5 for all turbulent simulations). 



4 Numerical experiments 

All numerical experimental compare the results of a GKS and an Navier-Stokes 
scheme. Both schemes are implemented into a 2D structured, finite-volume 
spatial discretization. The two schemes share the reconstruction of the conser- 
vative variables at cells interfaces, but differ in the evaluation of fluxes and in 
time stepping. Whereas the Navier-Stokes fluxes are obtained from Roe's ap- 
proximate Riemann solver (advective) and from central differences (viscous) , the 



GSK fluxes are obtained from Xu's scheme 19 extended to multi-dimensions. 
Navier-Stokes are advanced by means of a third order RK whereas GKS uses 
a time-accurate single-step approach. Both schemes use pre-conditioning (ap- 
proximate LU-SGS based on the approximate factorization of the Navier-Stokes 
operator [22] plus local time-stepping and multigrid acceleration ( |9j). The 
approximate LU-SGS factorization had already been used with GKS operator 



in 20 . Most computations have been conducted with CFL > 5. 

The evaluation of GKS fluxes requires roughly three times longer than the 
Navier-Stokes. However, the evaluation of fluxes is required only once in a 
time step, whereas the Navier-Stokes requires multiple evaluations (depending 
on the scheme). Broadly speaking, the explicit schemes have very close time 
performances whereas the use of pre-conditioning makes the Navier-Stokes ap- 
proximately twice faster. 
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4.1 RAE2822 and NACA 0012 airfoil in transonic flow 



Cases 9 and 10 of the measurements conducted by Cook ( [6]) on the RAE2822 
supercritical airfoil as well as Harris' investigation ( [8]) of the transonic flow 
around the NACA 0012 airfoil are arguably the most popular benchmarks for 
Navier-Stokes solvers and turbulence models, developed for transonic flow. In 
case 9 the boundary layer does not separate whereas in Case 10 and around the 
NACA 0012 at M = 0.800, Re = 9 x 10 6 and a = 2.83° the shock - boundary 
layer interaction leads to a large flow separation on the upper side of the airfoil. 

It is well known that two-equation turbulence models, such as the k — ui, 
become less and less accurate as the separated region grows. In the two sep- 
arated cases shown here, the size of the separated area is typically underesti- 
mated and the position of the shock predicted slightly downstream (refer for 
instance to [l7])- Figures [l] and [2] show the pressure distribution obtained with 
the Navier-Stokes and the GKS schemes. The computational meshes are C- 
type with 625 x 125 cells (RAE 2822) and 624 x 128 cells (NACA 0012), with 
Di<l (resolution of the boundary layer in wall units) in both cases. Results 
of Navier-Stokes and GKS are comparable for Case 9 (figure [I] lhs), whereas 
the GKS shows a higher accuracy in the capture of the shock - boundary layer 
interaction in Case 10 and in the NACA 0012 case (figure [I] rhs and figure [2] re- 
spectively), in that the flow separation is predicted much more accurately. The 
poor prediction of separated flows by the k — uj model can be improved by the 
replacing the linear expression for the eddy viscosity with algebraic relations for 
the components of the Reynolds stresses tensor ( [Tt] ) . 



4.2 Airfoil NACA 64A010 in transonic flow at high angle 
of attack 

The serie-6 airfoil NACA 64A010 has been investigated in transonic flow by 



Johnson [To]. The case Re = 2.0 x 10 6 , M = 0.75, angle of attack a = 6.2° has 
been considered here and investigated with the Navier-Stokes and GKS schemes. 
This flow case include a large separation with the shock wave located at about 
30% of chord. Neither the Navier-Stokes nor the GKS manage to capture the 
shock position accurately (figure [3]). However, this case provides the evidence 
of a different prediction of vortical / turbulent flow. The re-circulation area 
by the trailing edge is predicted to have two different patterns by the Navier- 
Stokes (figure [I| and GKS (figure [5]), especially around the trailing edge (figure 
[6]). Surprisingly, Johnson envisages re-circulation pattern similar to the one 
predicted by the GKS scheme ( [l0])- The computational meshes are C-type 
with 576 x 128 cells, and a resolution of the boundary layer of < 1 in wall 
units. 
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Gas-Kinetic Scheme 

Navier-Stokes (Roe) — - 
Experi menial Data (Cook 1981) 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

x/chord 



Gas-Kinetc Scheme - 
Navier-Stokes (Roe) - 
mental Data (Cook 1981) 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

x chord 



Figure 1: Airfoil RAE2822, Case 9 (lhs) Re = 6.2 x 10 6 , M = 0.725, angle of 
attack a = 2.30°. Case 10 (rhs) Re = 6.2 x 10 6 , M = 0.745, angle of attack 
a = 2.30°. Pressure coefficient computed and measured. 



5 Conclusions 

The GKS investigated in this study ( [19]) provides two main advantages with 
respect to ordinary Navier-Stokes schemes: higher-order fluxes and simultaneous 
treatment of transport and collisions. In the numerical experiments described, 
the GKS has in fact outperformed a Navier-Stokes scheme in a number of flow 
cases, characterized by interaction between shock and turbulent boundary layer. 
Interestingly, the resolution of secondary flows suggests that the GKS treats 
turbulence in a way similar to higher-order turbulence models (e.g. [17]), which 
are based on assumptions on the type of turbulence. 

The turbulent GKS seem to be a good candidate to investigate turbulent 
flow and the more so in the rarehed / transition regime, such as for instance 
flow cases related to hypersonic flight, where ordinary schemes still fail to pro- 
vide accurate and reliable results ( (|14j[T5]). Not only are GKS much better 
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Gas- Kinetic Scheme - 
Navier-Stokes (Roe) - 
al Data(Harris 1979) 



Figure 2: Airfoil NACA 0012, Re = 9.0 x 10 6 , M = 0.799, angle of attack 
a = 2.26°. Pressure coefficient computed and measured. 




Figure 3: Airfoil NACA 64A010, Re = 2.0 x 10 6 , M = 0.75, angle of attack 
a = 6.2°. Pressure coefficient computed and measured. 



Figure 4: Airfoil NACA 64A010, Re = 2.0 x 10 6 , M = 0.75, angle of attack a = 
6.2°. Streamlines showing the large separation induced by the shock obtained 
from the Navier-Stokes solutions. 
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Figure 5: Airfoil NACA 64A010, Re = 2.0 x 10 6 , M = 0.75, angle of attack a = 
6.2°. Streamlines showing the large separation induced by the shock obtained 
from the the Navier-Stokes (lhs) and GKS (rhs) solutions. Details of streamlines 
around the trailing edge. 




Figure 6: Airfoil NACA 64A010, Re = 2.0 x 10 6 , M = 0.75, angle of attack 
a = 6.2°. Streamlines showing the detail of the secondary flow around the 
trailing edge, obtained from the GKS solution. 

suited than Navier-Stokes schemes to handle flows with not negligible Kn, but 
they would also provide advantages in turbulence modelling. It is worth re- 
minding that virtually all turbulence theories have been developed under the 
assumption that the turbulent timescales are much bigger than molecular ones. 
As a matter of fact, in a flow characterized by M = 10 and Kn = 0.001 the 
two timescales might be comparable, i.e. r t ~ r. This means that ordinary 
Navier-Stokes schemes not only separate transport and collisions but they also 
miss the interactions between molecular and turbulent dynamics. 

Even in the continuum regime the properties of GKS schemes are much less 
known than those of Navier-Stokes schemes: many aspects of GKS still need 
to be clarified and future activities might include the sensitivity of results to 
the type of turbulence model, reconstruction order and truncation order in the 
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Chapman-Enskog expansion. 
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